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Abstract 

Determinantal Point Processes (Dpps) are elegant probabilistic models of repulsion and 
diversity over discrete sets of items. But their applicability to large sets is hindered by expensive 
cubic-complexity matrix operations for basic tasks such as sampling. In light of this, we propose 
a new method for approximate sampling from discrete k-D r i’s. Our method takes advantage 
of the diversity property of subsets sampled from a Dpp, and proceeds in two stages: first it 
constructs coresets for the ground set of items; thereafter, it efficiently samples subsets based on 
the constructed coresets. As opposed to previous approaches, our algorithm aims to minimize 
the total variation distance to the original distribution. Experiments on both synthetic and real 
datasets indicate that our sampling algorithm works efficiently on large data sets, and yields 
more accurate samples than previous approaches. 


i Introduction 

Subset selection problems lie at the heart of many applications where a small subset of items 
must be selected to represent a larger population. Typically, the selected subsets are expected 
to fulfill various criteria such as sparsity, grouping, or diversity. Our focus is on diversity, a 
criterion that plays a key role in a variety of applications, such as gene network subsampling [9], 
document summarization [36], video summarization [23], content driven search [4], recommender 
systems [47], sensor placement [29], among many others [1, 5, 21, 30, 33, 43, 44]. 

Diverse subset selection amounts to sampling from the set of all subsets of a ground set according 
to a measure that places more mass on subsets with qualitatively different items. An elegant 
realization of this idea is given by Determinantal Point Processes (Dpps), which are probabilistic 
models that capture diversity by assigning subset probabilities proportional to (sub)determinants 
of a kernel matrix. 

Dpps enjoy rising interest in machine learning [4, 23, 28, 30, 32, 34, 38]; a part of their appeal 
can be attributed to computational tractability of basic tasks such as computing partition functions, 
sampling, and extracting marginals [27, 33]. But despite being polynomial-time, these tasks remain 
infeasible for large data sets. Dpp sampling, for example, relies on an eigendecomposition of the 
Dpp kernel, whose cubic complexity is a huge impediment. Cubic preprocessing costs also impede 
wider use of the cardinality constrained variant /c-D p p [32]. 

These drawbacks have triggered work on approximate sampling methods. Much work has been 
devoted to approximately sample from a Dpp by first approximating its kernel via algorithms such 
as the Nystrom method [3], Random Kitchen Sinks [2, 41], or matrix ridge approximations [45, 
46], and then sampling based on this approximation. However, these methods are somewhat 
inappropriate for sampling because they aim to project the Dpp kernel onto a lower dimensional 
space while minimizing a matrix norm, rather than minimizing an error measure sensitive to 
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determinants. Alternative methods use a dual formulation [30], which however presupposes a 
decomposition L = XX T of the DPP kernel, which may be unavailable and inefficient to compute 
in practice. Finally, MCMC [6, 10, 14, 28] offers a potentially attractive avenue different from the 
above approaches that all rely on the same spectral technique. 

We pursue a yet different approach. While being similar to matrix approximation methods in 
exploiting redundancy in the data, in sharp contrast to methods that minimize matrix norms, we 
focus on minimizing the total variation distance between the original Dpp and our approximation. 
As a result, our approximation models the true Dpp probability distribution more faithfully, while 
permitting faster sampling. We make the following key contributions: 

- An algorithm that constructs coresets for approximating a k-D r p by exploiting latent struc¬ 
ture in the data. The construction, aimed at minimizing the total variation distance, takes 
0 (NM 3 ) time; linear in the number N of data points. The construction works as the overhead 
of sampling algorithm and is much faster than standard cubic-time overhead that exploits 
eigendecomposition of kernel matrices. We also investigate conditions under which such an 
approximation is good. 

- A sampling procedure that yields approximate /c-D r P subsets using the constructed coresets. 
While most other sampling methods sample diverse subsets in 0 (k 2 N) time, the sampling 
time for our coreset-based algorithm is 0 (k 2 M), where M C N is a user-specified parameter 
independent ofN. 

Our experiments indicate that our construction works well for a wide range of datasets, delivers 
more accurate approximations than the state-of-the-art, and is more efficient, especially when 
multiple samples are required. 

Overview of our approach. Our sampling procedure runs in two stages. Its first stage constructs 
an approximate probability distribution close in total variation distance to the true k-D pp. The 
next stage efficiently samples from this approximate distribution. 

Our approximation is motivated by the diversity sampling nature of Dpps: in a Dpp most of the 
probability mass will be assigned to diverse subsets. This leaves room for exploiting redundancy. 
In particular, if the data possesses latent grouping structure, certain subsets will be much more 
likely to be sampled than others. For instance, if the data are tightly clustered, then any sample 
that draws two points from the same cluster will be very unlikely. 

The key idea is to reduce the effective size of the ground set. We do this via the idea of 
coresets [17, 25], small subsets of the data that capture function values of interest almost as well as 
the full dataset. Here, the function of interest is a k-Dvv distribution. Once a coreset is constructed, 
we can sample a subset of core points, and then, based on this subset, sample a subset of the 
ground set. For a coreset of size M, our sampling time is 0 (k 2 M), which is independent ofN since 
we are using k-Dpps [32]. 

Related work. Dpps have been studied in statistical physics and probability [11, 12, 27]; they have 
witnessed rising interest in machine learning [21, 23,30,32-34,38, 44]. Cardinality-conditioned Dpp 
sampling is also referred to as "volume sampling", which has been used for matrix approximations 
[15, 16]. Several works address faster Dpp sampling via matrix approximations [3, 15, 30, 37] or 
MCMC [10, 28]. Except for MCMC, even if we exclude preprocessing, known sampling methods 
still require 0 (k 2 N) time for a single sample; we reduce this to 0 (k 2 M). Finally, different lines of 
work address learning DPPs [4, 22, 31, 38] and MAP estimation [20]. 

Coresets have been applied to large-scale clustering [8, 18, 24, 26], PCA and CCA [18, 39], and 
segmentation of streaming data [42]. 
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2 Setup and basic definitions 

A determinantal point process Dpp(L) is a distribution over all subsets of a ground set y of 
cardinality N. It is determined by a positive semidefinite kernel L G 1 R NxN . Let Ly be the submatrix 
of L consisting of the entries L,y with i,j G Y C y. Then, the probability P^(Y) of observing Y C y 
is proportional to det(Ly); consequently, Pi(Y) — det(Ly)/ det(L + I). Conditioning on sampling 
sets of fixed cardinality k, one obtains a k -D P P [32]: 

PL,k(y) ■■ = Pl(Y I \Y\ =k) 

-det(L y ) e ,(L)- 1 [|Y| =k} t 

where e k (L) is the k-th coefficient of the characteristic polynomial det(AI — L) = J^L 0 (—l) k e k (L)A N ~ k . 
We assume that Pl,1c(Y) > 0 for all subsets Y C y of cardinality k. To simplify notation, we also 
write P k = P Lrk . 

Our goal is to construct an approximation P k to P k that is close in total variation distance 

\\Pk-Pk\\tv:=l E \Pk(y)-Pk(Y)\, (2.1) 

Ycy, |y|=it 

and permits faster sampling than P k . Broadly, we proceed as follows. First, we define a partition 
TI — {Yi ,..., y M ) of y and extract a subset C C y of M core points, containing one point from 
each part. Then, for the set C we construct a special kernel L (as described in Section 3). When 
sampling, we first sample a set Y ~ Dpfj-(L) and then, for each c G Y we uniformly sample one 
of its assigned points y G y c . These second-stage points y form our final sample. We denote the 
resulting distribution by P k = Pq k . Algorithm 1 formalizes the sampling procedure, which, after 
one eigendecomposition of the small matrix L. 

We begin by analyzing the effect of the partition on the approximation error, and then devise an 
algorithm to approximately minimize the error. We empirically evaluate our approach in Section 6. 


3 Coreset sampling 

Let TI — {Yi,..., y m} be a partition of y, i.e., U-^Y; = y and Yt n Yy = 0 for i 7^ j. We call C C y 
a coreset with respect to a partition LI if \C n = 1 for i G [M]. With a slight abuse of notation, 
we index each part y c G TI by its core c G C Hi 34 - Based on the partition TI, we call a set Y C y 
singular 1 with respect to TL C n, if for y, G LI' we have Y fl Y/ < 1 and for Y, G n\n / we have 
|Y D Y/| = 0 . We say Y is k-singular if Y is singular and |Y| — k. 

Given a partition TI and core C, we construct a rescaled core kernel L G TR ' VI x M with entries 
L cc ' — \J | y : c 11 3 v I L c ■ We then use this smaller matrix L and its eigendecomposition as an input 
to our two-stage sampling procedure in Algorithm 1, which we refer to as CoreDpp. The two 
stages are: (i) sample a k-subset from C according to Dpp k (L); and (ii) for each c, pick an element 
y G y c uniformly at random. This algorithm uses only the much smaller matrix L and samples 
a subset from y in 0 (k 2 M) time. When M N and we want many samples, it translates into a 
notable improvement over the 0 (k 2 N) time of sampling directly from DpiyfL). 

The following lemma shows that CoreDpp is equivalent to sampling from a /c-Dpp where we 
replace each point in y by its corresponding core point, and sample with the resulting induced 
kernel Lc(y)- 

Lemma 1. CoreDpp is equivalent to sampling from Dpp k {L c ^), where in L C (yj we replace each 
element in y c by c,for all c G C. 

‘In combinatorial language, Y is an independent set in the partition matroid defined by It. 
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Algorithm 1: CoreDpf Sampling 

Input: core kernel L G and its eigendecomposition; partition FI; size k 

sample C Dpp t(L) 

sample y, ~ Uniform ( 34 ) for c G C 

return Y= {1/1,... ,y fc } 


Proof. We denote the distribution induced by Algo. 1 by Pc /k and that induced by D p by 

p’ 

1 k- 

First we claim that both sampling algorithms can only sample /c-singular subsets. By construc¬ 
tion, P C/k picks one or zero elements from each 34 - For Pj. r if Y is k-nonsingular, then there would 
be identical rows in (L C (jq)y = L C (y)> resulting in det(L C (y^) = 0 . Hence both Pc :k and P r k only 
assign nonzero probability to /c-singular sets Y. As a result, we have 

ek(Lc {y) )= E 011^1) det(Lc) 

C is k-singular cEC 

= E (II 1^1 det(Lc)) - E det ( £ c) - e k (L). 
ccc,|c|=fc cec |C|=Jfc 


For any Y = {\j\,.. ■ ,y k } C y that is k-singular, we have 


Pc,k(Y) 


det(L C ( Y )) 

ek(L)nU\y C ( yi )\ 


d et(k c(y) ) 

e k( L c(y)) 




(YlUiy^DdetjLc^) 
e k{Lc{y)) nf=l I yC (y/) I 


which shows that these two distributions are identical, i.e., sampling from Dr v k ( V) followed by 
uniform sampling is equivalent to directly sampling from Dpp^L^-yj). □ 


4 Partition, distortion and approximation error 

Let us provide some insight on quantities that affect the distance \\Pck ~ f\||tv when sampling 
with Algo. 1. In a nutshell, this distance depends on three key quantities (defined below): the 
probability of nonsingularity Su, the distortion factor 1 + ejy, and the normalization factor. 

For a partition n we define the nonsingularity probability S\ [ as the probability that a draw 
Y ~ Dppj-(L) is not singular with respect to any TL C n. 

Given a coreset C, we define the distortion factor 1 + £n (for £n > 0 ) as a partition-dependent 
quantity, so that for any c G C, for all u, v G 34 , and for any (k — l)-singular set S with respect to 
n \ 34 the following bound holds: 

de t(bsu{w}) Lu,u ~ L u $Lc Lg u ^ . . 

--77-—-r = - £ -— < 1 + £n- ( 4 -l) 

d et(U su{ ^ } ) L V ' V -L ViS Lf x L SiV ~ 

If f is the feature map corresponding to the kernel L, then geometrically, the numerator of (4.1) is 
the length of the projection of f (u) onto the orthogonal complement of span{ f(s) | s G S}. 

The normalization factor for a 4 -Dpp (L) is simply e k (L). 

Given n, C and the corresponding nonsingularity probability and distortion factors, we have 
the following bound: 
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Lemma 2. Let Y ~ Dppj-(L) and C{Y) be the set where we replace each y G Y by its core c G C, i.e., 
y G y c - With probability 1 - it holds that 


(l + en)- fc < 


det(L C ( Y )) 

det(Ly) 


< (l + £n)*. 


( 4 - 2 ) 


Proof. Let c G C and consider any (fc — l)-singular set S with respect to II \ y c . Then, for any 
v G y c , using Schur complements and by the definition of £n we see that 


- ^ de t(^su{c}) _ Lc,c Lg /C 

£n “ det(L suW ) “ Lv , v - L ViS L~ s l L SiV 


IIQ S ^( £ )II 2 

WQs^W 


< (1 + £ n)- 


Here, Q s _l is the projection onto the orthogonal complement of span{<^(s) | s G S}, and <p the 
feature map corresponding to the kernel L. 

With a minor abuse of notation, we denote by C{y) — c the core point corresponding to y, 
i.e., y G y c . For any Y = {y lr .. .,y k }, we then define the sets Y t = (C(j/i), ... ,C(yi),y i+1 ,.. .,y k }, 
where we gradually replace each point by its core point, with Yo — Y. If Y is /c-singular, then 
C (y;) C (yj) whenever i 7^ j, and, for any 0 < i < k — 1 , it holds that 


(1 + e n ) 1 < 


det(Ly. +1 ) 

det(Ly.) 


< 1 + £ n - 


Hence we have 


(1 + £ n) -fc 


< det(L C ( Y )) 

— det(Ly) 


rj det ( L y, +1 ) 

i = 0 det(Ly.) 


< (l + £ n) fc - 


This bound holds when Y is /c-singular, and, by definition of [, this happens with probability 

1 - < 5 n - □ 

Assuming £n is small. Lemma 2 states that if replacing a single element in a given subset with 
another one in the same part does not cause much distortion, then replacing all elements in the 
subset with their corresponding cores will cause little distortion. This observation is key to our 
approximation: if we can construct such a partition and coreset, we can safely replace all elements 
with core points and then approximately sample with little distortion. More precisely, we then 
obtain the following result that bounds the variational error. Our subsequent construction aims to 
minimize this bound. 


Theorem 3. Let P k — Dpp^.(L) and let be the distribution induced by Algo. 1. With the normalization 
factors Z = e k (L) and Zq = e k (L), the total variation distance between P k and Pq^ is bounded by 

\\Pk ~ Ffc^Htv < |1 — %-| + he n + (1 - k£n)^n- 

Proof. From the definition of Z and Zq we know that Z — 53 |y|=fc det(Ly) and 

Z C = det (( L C(T))Y) = X] det (Fc(y)) 

\Y\=k 

— J2 det (^C(Y))- 
Y ^-singular 
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The last equality follows since, as argued above, detftyyy)) = 0 for nonsingular Y. It follows that 

ll P fc- p C,Jcl|tv = E \Pk(V) - PcA Y )\ 

|Y| =k 

= E \Pk(r)-Pc*(Y)\+ E fltQO- ( 4 - 3 ) 

Y ^-singular Y A:-nonsingular 


For the first term, we have 


E \ p k(y) - PcA Y )\ = E 

Y /c-singular Y /c-singular 


det(Ly) det(L C ( Y )) 


< E 

Y /c-singular 

= \ E det(L Y ) 

Y /c-singular 

< ke u ( 1 - t> n ) + 


-(det(Ly) -det(L c(y) )) 


1 - 


det(L C (y)) 


det(Ly) 


+ E 

Y /c-singular 

i i 
z “ 


det(L C( yq) 


1 

Zc 


+ Z C 


’-I 


where the first inequality uses the triangle inequality and the second inequality relies on Lemma 2. 
For the second term in (4.3), we use that, by definition of 5 n, 

E Pk{Y) - *n- 

Y /c-nonsingular 

Thus the total variation difference is bounded as 


I Pt-P, 


C,k 11 tv 2s 


< 


'“I 

'“I 


+ fe n (l - <$n) +^n 
+ k£ U + (1 -fce n )(5n- 




In essence, if the probability of nonsingularity and the distortion factor are low, then it is 
possible to obtain a good coreset approximation. This holds, for example, if the data has intrinsic 
(grouping) structure. In the next subsection we provide further intuition on when we can achieve 
low error. 


4.1 Sufficient conditions for a good bound 

Theorem 3 depends on the data and the partition IT. Here, we aim to obtain some further intuition 
on the properties of n that govern the bound. At the same time, these properties suggest sufficient 
conditions for a "good" coreset C. For each y c , we define the diameter 

pc * = max \/L llu T L vv 2 L UZ7 . ( 4 - 4 ) 

u ,vey c 

Next, define the minimum distance of any point u G 34 to the subspace spanned by the feature 
vectors of points in a "complementary" set S that is singular with respect to n \ 34 > 

: = \l de det(L s f = \/ Lu ' u ~ L u,S L ~S lL S,u- 

Lemma 4 connects these quantities with £n; it essentially poses a separability condition on n (i.e., 
n needs to be "aligned" with the data) so that the bound on £n holds. 
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Lemma 4. If d c > p c for nil c € C, then 


£n < max 

ceC 


(2 d c 
(d c 


Pc)Pc 

Pc ) 2 


( 4 - 5 ) 


Proof. For any c £ C and any u, v & 34 and S (k — l)-singular with respect to IVfy c , we have 
det(L SU | u |) _ det(L s )(L lhU - L lhS Lg 1 L S/U ) 

det(Ls)(L VfV L U; sLg Lg /V .) 

_ Lu,u — L u sLg Lgu _ IIQs ± ^ , ( u ) IP 

~ L V/V - L^L^L^, _ JQ^W 

Without loss of generality we assume det(L su r,q) > det(L Sj j l ,|). By definition of p c we know that 


0 < HQ s x^(m)|| - ||Q s ±^(p)|| 

< ||Q s x (</>(«)-</>(»)) II < \\<P( U ) ~<P( V )\\ — Pc * 

Since 0 < ||Q s ± 0 (p)|| < ||Q s _l$(u)|| < ||<p(u)|| by assumption, we have 

llQs^OOII 2 IIQ S ^(»)|| 2 

||Q S ±<K^)II 2 “ (IIQs-l^(w)II ~pc ) 2 

< ( IlftOOll \ 2 < ( d c \ 2 
^'\\\(p(u)\\-p c J ~ \(d c -p c )J ‘ 


Then, by definition of £n, we have 


1 + £n < max 

C 


d 2 

(d c ~ pc) 2 ' 


from which it follows that 


£ n < max 

C 


(2 d c pc)pc 
(d c -pc) 2 


□ 


5 Efficient construction 

Thm. 3 states an upper bound on the error induced by CoreDpp and relates the total variation 
distance to n and C. Next, we explore how to efficiently construct n and C that approximately 
minimize the upper bound. 

5.1 Constructing TT 

Any set Y sampled via CoreDpp is, by construction, singular with respect to n. In other words, 
CoreDpp assigns zero mass to any nonsingular set. Hence, we wish to construct a partition n 
such that its nonsingular sets have low probability under Dpp^(L). The optimal such partition 
minimizes the probability 6 [ \ of nonsingularity. A small < 5 n value also means that the parts of n 
are dense and compact, i.e., the diameter p c in Equation (4.4) is small. 

Finding such a partition optimally is hard, so we resort to local search. Starting with a current 
partition n, we re-assign each y to a part y c to minimize < 5 n- If we assign y to y c , then the 
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probability of sampling a set Y that is singular with respect to the new partition n is 

\ 

P[Y ~ DpP jt(L) is singular] — — ^ det(Ly) 

Y/c-singular 

= 4( L det(Ly)+ £ det (M) 

Yk-sing.,y£Y Yk-sing.,yeY 

= 4 (const + d et(b ru{y} )) 

Y' (k- l)-sing. w.r.t n \ y c 

= \ ( const + L yy s k-i ( L l)) / 

where s^(L) 5 Zy fc-sing. det(Ly). The matrix L\ c denotes L with rows y c and columns 34 deleted, 
and U — L — Ly^L^y. For local search, we would hence compute Li /1/ s[ I _ | (jX c ) for each point y 
and core c, assign y to the highest-scoring c. Since this testing is still expensive, we introduce 
further speedups in Section 5.3. 

5.2 Constructing C 

When constructing C, we aim to minimize the upper bound on the total variation distance between 
Pj c and Pq k stated in Theorem 3. Since < 5 n and e rt only depend on n and not on C, we here focus 
on minimizing |1 — |, i.e., bringing Zq as close to Z as possible. To do so, we again employ local 

search and subsequently swap each c € C with its best replacement v & y c . Let C c,v be C with c 
replaced by v. We aim to find the best swap 

v = argmin t , e3 ; c |Z - Z \ ( 5 - 1 ) 

= argmin^g-yJZ - e k (L C c,v (y) )\. (5.2) 

Computing Z requires computing the coefficients tyf Lj, which takes a total of 0 (N 3 ) time 2 . In the 
next section, we therefore consider a fast approximation. 

5.3 Faster constructions and further speedups 

Local search procedures for optimizing LI and C can be further accelerated by a sequence of 
relaxations that we found to work well in practice (see Section 6). We begin with the quantity 
sF_i(Lf c ) that involves summing over sub-determinants of the large matrix L. Assuming the 
initialization is not too bad, we can use the current C to approximate y. In particular, when 
re-assigning y, we substitute all other elements with their corresponding cores, resulting in the 
kernel L — ^ciyj ■ This changes our objective to finding the c £ C that maximizes s^ 1 (Lf c ). Key to 
a fast approximation is now Lemma 5, which follows from Lemma 1. 

Lemma 5. For all k < |TT|, it holds that 

s V( L C(y)) = e k(L C (y)) = e k (L). 

2 In theory, this can be computed in 0(N U> log (N)) time [ 13 ], but the eigendecompositions and dynamic programming 

used in practice typically take cubic time. 
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Algorithm 2: Iterative construction of n and C 

Require: II initial partition; C initial coreset; k the size of sampled subset; v number of nearest neighbors 
taken into consideration 
while not converged do 

for all y 6 y do 

c <— group in which y lies currently: y 6 34 

if y € C then 

continue 

end if 

G {groups of v cores nearest to X v } 
g* = argmax^gS^OL^) 

if c / g* then 

34 = 34 \{y} 

y? = 3V u {y} 

if e fc( L cs*hy)) > e ki L c'>i(y)) then 

c <- cf* 

end if 
end if 
end for 

for all g 6 [M] do 

; = arg maXj e y g e k{L cg .j(y)) 

C = CS’i 

end for 
end while 


Proof. 

s ¥(^c(y)) = det ((^c(y))Y) — d e t(L C ( Y )) 

Y k- sing. Y k- sing. 

— det ( L c(y)) = e k( L c(y)) — e k(Q; 

\Y\=k 

the last equality was shown in the proof of Thm. 3. □ 

Computing the normalizer ejt(L) only needs 0 (M 3 ) time. We refer to this acceleration as 
CoreDpp-z. 

Second, when constructing C, we observed that Zq is commonly much smaller than Z. Hence, 
a fast approximation merely greedily increases Zq without computing Z. 

Third, we can be lazy in a number of updates: for example, we only consider changing cores 
for the part that changes. When a part y c receives a new member, we check whether to switch 
the current core c to the new member. This reduction keeps the core adjustment at time 0 (M 3 ). 
Moreover, when re-assigning an element y to a different part 34 , it is usually sufficient to only 
check a few, say, v parts with cores closest to y, and not all parts. The resulting time complexity for 
each element is 0 (M 3 ). 

With this collection of speedups, the approximate construction of n and C takes 0 (NM 3 ) for 
each iteration, which is linear in N, and hence a huge speedup over direct methods that require 
0 (N 3 ) preprocessing. The iterative algorithm is shown in Algorithm 2. The initialization also 
affects the algorithm performance, and in practice we find that kmeans++ as an initialization works 
well. Thus we use CoreDpp to refer to the algorithm that is initialized with kmeans++ and 
uses all the above accelerations. In practice, the algorithm converges very quickly, and most of 
the progress occurs in the first pass through the data. Hence, if desired, one can even use early 
stopping. 
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Figure 1 : Total variation distances (error) on synthetic data with varying nChist and /'2-norm. 


6 Experiments 

We next evaluate CoreDpp, and compare its efficiency and effectiveness against three competing 
approaches: 

- Partitioning using /c-means (with kmeans++ initialization [7]), with C chosen as the centers of 
the clusters; referred to as K++ in the results. 

- The adaptive, stochastic Nystrom sampler of [3] ( NysStoch ). We used M dimensions for NysStoch, 
to use the same dimensionality as CoreDpp. 

- The Metropolis-Hastings DPP sampler MCDPP [28]. We use the well-known Gelman and Rubin 
multiple sequence diagnostic [19] to empirically judge mixing. 

In addition, we show results using different variants of CoreDpp: CoreDpp-z described 
in Sec. 5.3 and variants that are initialized either randomly (CoreDpp-r) or via kmeans++ (CoreDpp). 

6.1 Synthetic Dataset 

We first explore the effect of our fast approximate sampling on controllable synthetic data. The 
experiments here compare the accuracy of the faster CoreDpp from Section 5.3 to CoreDpp-z, 
CoreDpp-r and K++. 

We generate an equal number of samples from each of nClust 30-dimensional Gaussians with 
means of varying length (/2-norm) and unit variance, and then rescale the samples to have the 
same length. As the length of the samples increases, £n and d\\ shrink. Finally, L is a linear kernel. 
Throughout this experiment we set k = 4 and N — 60 to be able to exactly compute II Pk ~ P/clltv- 
We extract M = 10 core points and use v — 3 neighboring cores. Recall from Sec. 5.3 that when 
considering the parts that one element should be assigned to, it is usually sufficient to only check v 
parts with cores closest to y. Thus, v — 3 means we only consider re-assigning each element to its 
three closest parts. 

Results. Fig. 1 shows the total variation distance ||Pj- -Pfclltv defined in Equation (2.1) for the 
partition and cores generated by K++, CoreDpp, CoreDpp-r and CoreDpp-z as nClust and the 
length vary. We see that in general, most approximations improve as £n and />n shrink. Remarkably, 
the CoreDpp variants achieve much lower error than K++. Moreover, the results suggest that the 
relaxations from Section 5.3 do not noticeably increase the error in practice. Also, CoreDpp-r 
performs comparable with CoreDpp, indicating that our algorithm is robust against initialization. 
Since, in addition, the CoreDpp construction makes most progress in the first pass through 
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Figure 2: Approximate total variation distances (empirical estimate) on MNIST (left) and 
GENES (right) with M varying from 20 to 100 and fixed k — 6 . 




k k 

Figure 3: Approximate total variation distances (empirical estimate) on MNIST (left) and 
GENES (right) with k varying from 8 to 2 and fixed M = 100 . 


the data, and the kmeans++ initialization yields the best performance, we use only one pass of 
CoreDpp initialized with kmeans++ in the subsequent experiments. 

6.2 Real Data 

We apply CoreDpp to two larger real data sets: 

1. MNIST [35]. MNIST consists of images of hand-written digits, each of dimension 28 x 28 . 

2. GENES [9]. This dataset consists of different genes. Each sample in GENES corresponds to a 
gene, and the features are shortest path distances to 330 different hubs in the BioGRID gene 
interaction network. 

For our first set of experiments on both datasets, we use a subset of 2000 data points and an RBF 
kernel to construct L. To evaluate the effect of model parameters on performance, we vary M from 
20 to 100 and k from 2 to 8 and fix v — 2 (see Section 6.1 for an explanation of the parameters). 
Larger-scale experiments on these datasets are reported in Section 6.3. 
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Figure 4: (a) Example coreset C of size 40, each figure is a core in the coreset constructed by our 
algorithm; (b,c) Two different parts corresponding to the first and second core. 


Performance Measure and Results. On these larger data sets, it becomes impossible to com¬ 
pute the total variation distance exactly. We therefore approximate it by uniform sampling and 
computing an empirical estimate. 

The results in Figure 2 and Figure 3 indicate that the approximations improve as the number of 
parts M increases and k decreases. This is because increasing M increases the models' approxi¬ 
mation power, and decreasing k leads to a simpler target probability distribution to approximate. 
In general, CoreDpp always achieves lower error than K++, and NysStoch performs poorly in 
terms of total variation distance to the original distribution. This phenomenon is perhaps not so 
surprising when recalling that the Nystrom approximation minimizes a different type of error, a 
distance between the kernel matrices. These observations suggest to be careful when using matrix 
approximations to approximate L. 

For an intuitive illustration. Figure 4 shows a core C constructed by CoreDpp, and the elements 
of one part y c . 

6.3 Running Time on Large Datasets 

Lastly, we address running times for CoreDpp, NysStoch and the Markov chain /c-DPP (MCDPP [28]) 
For the latter, we evaluate convergence via the Gelman and Rubin multiple sequence diagnostic [19]; 
we run 10 chains simultaneously and use the CODA [40] package to calculate the potential scale 
reduction factor (PSRF), and set the number of iterations to the point when PSRF drops below 1.1. 
Finally we run MCDPP again for this specific number of iterations. 

For overhead time, i.e., time to set up the sampler that is spent once in the beginning, we 
compare against NysStoch: CoreDpp constructs the partition and L, while NysStoch selects 
landmarks and constructs an approximation to the data. For sampling time, we compare against 
both NysStoch and MCDPP: CoreDpp uses Algo. 1, and NysStoch uses the dual form of k-Dvv 
sampling [30]. We did not include the time for convergence diagnostics into the running time of 
MCDPP, giving it an advantage in terms of running time. 

Overhead. Fig. 5 shows the overhead times as a function of N. For MNIST we vary N from 6,000 
to 20,000 and for GENES we vary N from 6,000 to 10,000. These values of N are already quite 
large, given that the Dpp kernel is a dense RBF kernel matrix; this leads to increased running 
time for all compared methods. The construction time for NysStoch and CoreDpp is comparable 
for small-sized data, but NysStoch quickly becomes less competitive as the data gets larger. The 
construction time for CoreDpp is linear in N, with a mild slope. If multiple samples are sought, 
this construction can be performed offline as preprocessing as it is needed only once. 
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Figure 5: Overhead (setup) time in seconds with varying ground set size ( N ) on MNIST (left) and 
GENES (right). 




Figure 6: Average time for drawing one sample as the ground set size ( N ) varies on MNIST (left) 
and GENES (right). Note that the time axis is shown in log scale. 


Sampling. Fig. 6 shows the time to draw one sample as a function of N, comparing CoreDpp 
against NysStoch and MCDPP. CoreDpp yields samples in time independent of N and is extremely 
efficient - it is orders of magnitude faster than NysStoch and MCDPP. 

We also consider the time taken to sample a large number of subsets, and compare against both 
NysStoch and MCDPP —the sampling times for drawing approximately independent samples with 
MCDPP add up. Fig. 7 shows the results. As more samples are required, CoreDpp becomes 
increasingly efficient relative to the other methods. 


7 Conclusion 

In this paper, we proposed a fast, two-stage sampling method for sampling diverse subsets with 
k-D its. As opposed to other approaches, our algorithm directly aims at minimizing the total 
variation distance between the approximate and original probability distributions. Our experiments 
demonstrate the effectiveness and efficiency of our approach: not only does our construction have 
lower error in total variation distance compared with other methods, it also produces these more 
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Figure 7: Average time for sampling different numbers of subsets with N — 5000 , M — 40 and 
k = 5 on MNIST (left) and GENES (right). 


accurate samples efficiently, at comparable or faster speed than other methods. 
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